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Abstract 

We consider the numerical simulation of Hamiltonian systems of ordinary differen- 
tial equations. Two features of Hamiltonian systems are that energy is conserved 
along trajectories and phase space volume is preserved by the flow. We want to 
determine if there are integration schemes that preserve these two properties for 
all Hamiltonian systems, or at least for all systems in a wide class. This paper 
provides provides a negative result in the case of two dimensional (one degree of 
freedom) Hamiltonian systems, for which phase space volume is identical to area. 
Our main theorem shows that there are no computationally reasonable numerical 
integrators for which all Hamiltonian systems of one degree of freedom can be inte- 
grated while conserving both area and energy. Before proving this result we define 
what we mean by a computationally reasonable integrator. We then consider what 
obstructions this result places on the existence of volume- and energy-conserving 
integrators for Hamiltonian systems with an arbitrary number of degrees of freedom. 
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1 Introduction 



We consider a system of Hamiltonian differential equations on M. 2n 

dq dH dp dH 
dt dp ' dt dq ' 
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defined by the Hamiltonian function H : R 2ri — > R. We denote the time t flow 
map of these equations by S l . The flow of these differential equations has two 
important features. The first is that H is conserved along trajectories. That 
is, H(S\q,p)) = H(q,p), for all t G R and all (q,p) G R 2 ™. The second is that 
phase space volume is conserved by the flow: if A C R 2n is a bounded open 
set, then vol(S t A) = vo\(A) for all t. The latter property is a consequence of 
the symplecticity of the flow. (See, for example, [2]). 

For certain molecular dynamics applications, the ideal numerical integrator 
would retain these two properties of the flow [5] . That is, if the integrator with 
time step At defines a map $a< on R 2 ™, we would like both f/($At(?,£>)) = 
H(q,p) for all (q,p) (energy conservation) and vol(<E»At(^4)) = vol(A) for all 
bounded open subsets A of R 2ri (volume conservation). 

Symplectic integrators such as the implicit midpoint rule conserve volume ex- 
actly for all Hamiltonian systems, with any number of degrees of freedom, but 
do not conserve energy [2]. It has already been shown that in certain circum- 
stances it is unreasonable to expect that symplectic integrators which also 
conserve energy exist [6]. However, volume-conservation is a weaker property 
than symplecticity for Hamiltonian systems of more than one degree of free- 
dom (n > 1). Thus it seems plausible that there is a consistent integration 
scheme which for any Hamiltonian function H and time step At yields a map 
&At which conserves both volume and energy. 

In this article we will argue that there is no such integration scheme. We do this 
by showing that no numerical integrator is able to integrate all Hamiltonian 
systems in R 2 while simultaneously conserving energy and phase space volume. 
For special Hamiltonian systems in R 2 such an integrator is possible. (See [3] 
for examples of such systems of any number of degrees of freedom.) However, 
our theorem states that for any energy- conserving numerical integrator from 
a very broad class, there will be at least one Hamiltonian system on R 2 (and 
in fact very many) for which it does not conserve phase space volume. 

Before stating and proving our main theorem, in Section 2 we will define 
what we mean by an integrator. From the class of integrators we then define 
the computationally reasonable integrators. This will include any explicit or 
implicit formula for defining a new state (q k+1 ,p k+1 ) as a function of a state 
(q k ,p k ) and a timestep At such that the only information used about H is its 
value and the value of its derivatives at a finite number of points (r J , s J ) G 
R 2ra , j = 1, . . . , to. In Section 3 we will state and prove our main result in 
terms of this definition. We imagine that a numerical analyst has devised a 
computationally reasonable numerical integrator that is energy-conserving for 
all Hamiltonian systems in R 2 . We apply this integrator to the system with 
Hamiltonian function H(q, p) = p 2 /2 in R 2 . For two different inputs we observe 
at which points (r J , s- 7 ), j — 1, . . . , to the integrator depends on the function H 
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and its derivatives. Using this information, we construct another Hamiltonian 
function H(q,p) = p 2 /2 + V(q) which is arbitrarily close to H and has the 
following property: the numerical integrator cannot be volume-conserving for 
both H and H. 

The main result demonstrates that there is no integration scheme that con- 
serves volume and energy for arbitrary Hamiltonian systems of any number 
of degrees of freedom. However, we cannot conclude that there is no scheme 
that is energy- and volume-conserving for Hamiltonian systems in a particular 
number of dimensions 2n, n > 1. This question is open. To partially address 
this issue, in Section 4 we will show that under some reasonable — but not 
essential — conditions on an integrator, the problem reduces to that of the 
n = 1 case. Thus, there is no energy- and volume-conserving integrator for 
fixed n that satisfies these additional assumptions. 

We now explain the relation between the results here and those of the well- 
known paper of Zhong and Marsden [6]. In that paper the authors consider a 
Hamiltonian system for which there are no invariants except energy. They show 
that any integrator that is both symplectic and energy-conserving actually 
computes exact trajectories of the system up to a time reparametrization. 
From this result they conclude that energy-conserving symplectic integration 
is not possible in general, since presumably the set of Hamiltonians for which 
one could compute trajectories exactly up to a time reparametrization are very 
small. Though this argument is very plausible, it leaves open two questions: 

(1) How do we make precise the idea of something not being possible for any 
numerical integrator? 

(2) Is it possible to perform energy-conserving and symplectic integration for 
general Hamiltonians when we restrict ourselves to the IR 2 case? 

The importance of the latter question is that if it were possible, then volume- 
and energy- conserving integration would be possible for general Hamiltonian 
systems via J-splitting [3,1]. 

We provide an answer to the first question in Section 2 with the definition 
of a computationally reasonable integrator. As for the second question, since 
volume-conservation and symplecticity are identical in R 2 , Zhong and Mars- 
den's result shows that in R 2 volume- and energy-conserving integration is 
equivalent to solving the original system exactly up to a time-reparametrization 
(This result is stated and proved for this case as Lemma 3 in Section 3.) 
The main result of our paper answers the second question in the negative by 
showing that it is not possible for general Hamiltonian systems in R 2 using 
computationally reasonable integrators. 

Before we begin we discuss some interesting related work. Even though energy 
and volume conserving integrators may not exist, the paper [1] does the next 
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best thing. There the authors show how to approximate any Hamiltonian 
function arbitrarily well by a special piece-wise smooth function of a form 
described by [3] whose trajectories can be integrated while conserving volume 
and energy. The original Hamiltonian function is not conserved. However, 
unlike for standard symplectic methods, a modified Hamiltonian function close 
to the original is conserved exactly [1] for all time. 



2 What is a numerical integrator? 



An integrator for a Hamiltonian system of ordinary differential equations takes 
a Hamiltonian H, a step length At, and an initial value (q k ,p k ) G M 2n , and 
produces a value (q k+1 ,p k+1 ) G M. 2n . Typically, (q k+1 ,p k+1 ) depends on H 
through components of VH and perhaps H itself. If (q k ,p k ) is an approxi- 
mation to (q(kAt),p(kAt)), we take (q k+1 ,p k+1 ) to be an approximation to 
(q((k + l)At),p((k + l)At)). 

Definition 1 An integrator $ is a function that takes arguments H : R 2n — + 
R, At G [0, oo), (q k ,p k ) G R 2n and either returns (q k+1 ,p k+1 ) G R 2n or is not 
defined. We write 

(q k+ \p k+1 )=<S>(H, (q\p k ),At). 



We have allowed the integrator to not be defined for certain input values. This 
is often the case for implicit integrators when the vector field is insufficiently 
smooth or the time step is too large. 

In order to get meaningful constraints on what is computationally reasonable, 
we cannot let arbitrary maps be included in the class of algorithms we study. 
After all, the exact flow map has all the qualitative features one could want, 
but it is not feasible to compute it (even to machine precision) for most ap- 
plications. We would like our definition to be broad enough to include most 
existing numerical integrators. 

Informally, we say an integrator $ is computationally reasonable if for each 
(q k ,p k ) and At, (q k+1 ,p k+1 ) = (q k ,p k ),At) depends on H only through 

its value and the value of its derivatives at a finite number of points (r 7 , s- 5 ) G 
R 2n , j = 1, ... ,n. In the following formal definition we use multi-index no- 
tation to define higher-order derivatives: for H : R 2n — > R and a G Nq™ we 
let 

9aH := OqY ' ' ' dq^ dp" n+1 ' ' ' dp~^ H 
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Definition 2 An integrator <£> is computationally reasonable if for each H : 
R 2 ™ -> R, (q k ,p k ) G R 2ri and At > there exists 

(1) m e N 0; 

(2) (r*,si)eR 2n ,j = l,...,m, 

(3) a J GN^ = l,..,m, 

such that for any function H : R 2 ™ — > R £/ia£ satisfies 
d a jH(r j ) = d aj H(r j ), j = l,...,m 

eit/ier (q k ,p k ),At) = (q k ,p k ),At) or both are not defined. 

In the remainder of this section we discuss examples of the class of computa- 
tionally reasonable integrators. 

First note that all explicit methods fit into this class. We informally define 
an integrator to be explicit if it can be implemented by an algorithm that 
terminates in a finite number of steps using function evaluations of H or its 
derivatives, arithmetic operations, and logical operations. This includes all of 
the explicit Runge-Kutta and partitioned Runge-Kutta methods, for example. 
It also includes integrators that collect information adaptively to perform a 
step, such as the Bulirsch-Stoer method [4]. Taylor series methods are also 
included. 

The class of computationally reasonable integrators also includes implicit 
methods, such as the implicit Runge-Kutta methods. Here we define an im- 
plicit method to be one where (q k+1 ,p k+1 ) is specified by requiring it to be 
the solution to a nonlinear system of equations in H and its derivatives. Im- 
plicit algorithms cannot in general be implemented exactly in a finite number 
of steps, but they still fit into the framework of Definition 2. To see this, 
note that even though solving a system of nonlinear equations exactly typi- 
cally requires looking at H and its derivatives at an infinite number of points 
(while performing the Newton Iteration, for example), determining if we have 
a solution to a nonlinear system of equations only requires examining a finite 
number of points. So ifjg fc+1 ,p fc+1 ) solves the equations for H, it will still solve 
the equation for any H which is identical to H at the points. Similarly, step- 
and-project methods are included in this class. These are methods that consist 
of one step of a simpler method followed by a projection onto a manifold [2, 
IV.4]. 

What integrators do not satisfy Definition 2? Integrators that require the exact 
computation of integrals of H or its derivatives do not. Computing the integral 
of a general function requires knowing its value at an infinite set of points on 
the domain of integration. Unlike in the case of solving nonlinear equations, 
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if we are given a value for the integral, there is no finite set of points at which 
we can examine the function to verify that the value is correct. Of course, it 
is possible to define an numerical integrator that uses integrals of functions, 
but the actual computation of these integrals for general Hamiltonian systems 
would require numerical quadrature. This in turn would require sampling the 
function at a finite number of points, and introducing truncation error. The 
new method with this additional truncation error does form a computationally 
reasonable integrator, while the original method with the exact integral does 
not. 

Finally, we note that multistep methods are not even integrators according 
to Definition 1. We believe our framework could be extended to multistep 
method but we do not do so here. 



3 Main Theorem 

To prove our main result Theorem 5 we use the following lemma. It shows 
that for a Hamiltonian system in R 2 , the map defined by an energy and area 
preserving integrator is just a time-reparametrization of the flow map. As 
discussed in the introduction, this is essentially Zhong and Marsden's result 
[6] in the two dimensional case. The only addition is that we show that the 
time-reparametrization is just a constant rescaling of time where locally the 
constant does not depend on energy. 

Lemma 3 Let H : R 2 — > R be a smooth function and S f its induced Hamilto- 
nian flow map. Let E el be a particular energy. Let U be an open set whose 
intersection with {(q,p)\H(q,p) = E} is a simple curve T. Suppose that VH ^ 
on T. Let $At ; At > be a continuous area- conserving map defined on U 
that conserves H. Then there is a constant c such that $ At = S cAt (q,p) 
for all (q,p) € T. 

Proof: On U we can define canonical action-angle coordinates (9, 4>) in which 
the Hamiltonian function is H(9, <fi) = H{9). The flow map is then S t (9, 0) = 
(9, (f) + tH'(9)), where H'(9) ^ 0, since we still have [H'(ff), 0] T = VH ^ in 
the new coordinates. To study the Jacobian of $At i n the coordinates (9, 4>) 
we let (6>i,0i) = $a<(#o,0o) for (#o,0o) £ U. Then the Jacobian of <3>a< is 

1 

d<t>i dji 
_ 80 d<p . 

Area conservation implies that the determinant is 1 and so <90i/<90o — 1 for 
all 9q, 0o- This yields <fii = (fio + T for all 0o fc> r some r which may depend on 
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6q. Hence, if we choose c so that cAt/H'(6 ) = r we have 
$ At(0o, 0o) = (Oo, 0o + t) = S cAt (9 , O ), 

for all (#o, 0o) £ r, as required. O 

For the main result we will use the following very weak definition of consis- 
tency. 

Definition 4 Let S* be the flow map for the Hamiltonian H : R 2n — > R. An 

integrator $ is consistent at (q, p) £ R 2n if 

Um ||g A '(g,p)-$(ff,At,(g, P ))|| = Q 

At^O At 

Theorem 5 Let $ fre an integrator for a Hamiltonian system in R 2 . Let 
H(q,p) = p 2 /2. Let U be an open neighbourhood of {(q, l)\q £ R}. Suppose: 

(1) $ consistent at (0,1) (Definition 4)- 

(2) $ computationally reasonable (Definition 2). 

(3) $ conserves energy for any H , (q,p), and At for which it is defined. 

(4) For At < Ato, &{H, (q,p),At) is defined, depends continuously on (q,p), 
and conserves volume on U. 

Then for all sufficiently small At there is a function V : R — > R such 
that if H(q, p) = p 2 /2 + V(q) then ®(H, •, At) is not simultaneously defined, 
continuous, and volume-conserving on U. For each such At, the V constructed 
can be replaced by XV for A £ (0, 1) and the same result holds. 

Proof: For any At £ (0, Ato], since <$>(H, ■, At) conserves energy and is 
continuous on U, the set {(q, l)\q £ R} is mapped onto itself. As &(H, •, At) 
conserves volume, Lemma 3 shows that it is identical to the flow of the original 
Hamiltonian system on {(q, l)\q £ R} with a rescaling of time: 

*(H, (q,l),At) = (q + cAt,l), 

for all q £ R, where c does not depend q. The consistency condition at (0, 1) 
implies that for small enough At we have that c > 0. From now on, we assume 
At is small enough so that c > 0. 

Consider the integrator applied to H at the point (q,p) = (0,1). Since $ is 
computationally reasonable (Definition 2), there are a finite number of points 
(r J , s J ) £ R 2 , j = 1, . . . ,m, such that $ only depends on H at these points. 
Choose a q £ R big enough so that the interval [q , q + cAt] contains none 
of the points r J , j — 1, . . . , m and is disjoint from the interval [0, cAt]. 
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Consider the integrator applied to H at the point (go,l)- There are points 
(f- 7 , = 1, • • • , m such that &(H, (g , 1), At) only depends on H at these 
points. Let V be a function such that 

(1) V(g) = for q not in [g , g + cAt], and V(f j ) = 0, j = 1, . . . , m, 

(2) < V(q) < 1/2 for all g, 

(3) for some e > 0, V(g) > for g £ (g , go + e )- 

Note that multiplying V by any factor in (0, 1) gives a function satisfying the 
same conditions. 

Let H(q,p) = p 2 /2 + V(q) andjet S l denote the flow map of the system with 
this Hamiltonian. Now if &(H, -,At) were defined, continuous and volume- 
conserving on U, then by Lemma 3 

$(H, (q,l),At) = S dAt (q,l) 

for all q £ R for some constant d. We will show that this is impossible. 

First note that the integrator gives the same result for both Hamiltonians at 
(0, 1). This is because here the result of the integrator only depends on H only 
at the points (r 7 , s- 5 ), at which H and H agree. Since 

(dAt, 1) = 0>(H, (0, 1), At) = 0>(H, (0, 1), At) = (cAt, 1) 
we must have c = d. 

On the other hand, the integrators also give the same results for both Hamil- 
tonians at (go, 1)- This implies dm. the following way. We have that 

(g + cAt,l) = $(iJ,(g ,l),At) 
= $(H, (g ,l),At) 

= S dAt (q ,l) 
^S dAt (q ,l) 
= (q + dAt,l). 

The inequality follows from Lemma 6, since the Hamiltonians differ on the 
interval [go, go + e]. So c / d This is a contradiction. Therefore, •, At) 

cannot be simultaneously defined, continuous, and volume-conserving on U . 

O 

The following lemma asserts the intuitively clear fact that if V(q) is positive 
for q £ [go, go + e] then the flows starting from (g , 1) corresponding to H and 
H are different for positive times. 
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Lemma 6 Let V : R 2n ^ R be C°° with V(q ) = 0, 1/2 > V(q) > for 
all q G M and V(q) > /or all q G [go? <?o + e ]- H(q,p) = p 2 /2 and 
H(q,p) = p 2 /2 + V{q). Let S l and S l be the respective Hamiltonian flow maps 
ofH and H. Then S\q , 1) ^ 5*(g , 1) for t > 0. 

Proof. The trajectory for the Hamiltonian H as a function of time is (q(t),p(t)) = 
(t, 1) for all t. Letting q(t) describe the position for the Hamiltonian H, the 
usual solution technique gives 

f 1 
t = / dx =: F(q) 

I {l - 2V(x) 

The function F is strictly increasing and so has a well defined inverse. We can 
write q(t) = Now F'(q) > 1 for q > 0, so g'(i) < 1 for t > 0. Hence 

g(t) < t for t > and the two flow maps cannot be equal for t > 0. O 



4 Multiple Degrees of Freedom 



In the previous section we showed that there can be no general energy- and 
volume-conserving integration schemes because there are no integrators that 
conserve energy and volume for all Hamiltonians in M 2 . However, suppose we 
ask if such integrators exists for Hamiltonian systems of dimension 2n, n > 1. 
We conjecture that a result like Theorem 5 still holds in this case. However, 
the method of proof for the n — 1 case does not extend to this case. 

Instead we will state two conditions on an integrator, either one of which 
prevents it from being volume- and energy-conserving for general Hamiltonian 
systems in IR 2 ™. Both of these conditions are desirable for an integrator to have, 
but unlike computational reasonibility, it is not difficult to imagine a practical 
integrator that did not satisfy them. The proof of the theorems in this section 
will work by showing that, if a computational reasonable energy-conserving 
integrator with either condition exists, a special Hamiltonian in M 2n can be 
constructed for which it does not conserve volume. 

In the following, fix n > 2. Let q,p G lR n have components qi^pi G K. Define 
the projections 7Tj, % = 1, . . . , n by iTi(q,p) = (qi,Pi)- 

The first condition asserts that if variables (qi,Pi) do not occur in the Hamil- 
tonian function then the integrator does not change their value. 

Condition 7 If H(q,p) is independent of the variables (qi,pi) then 
n^(H, (q,p),At) = (qi,pi) 
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for all (q,p) and At for which <3> is defined. 

Theorem 8 Let $ be an integrator for Hamiltonian systems in M, 2n . Let 
H(q,p) — Pi/2. Let U be an open neighbourhood of the set 

{(( gi ,o,...,0),(l,0,...,0))| gi gK}. 
Suppose 

(1) $ is consistent (Definition 3). 

(2) $ is computationally reasonable (Definition 2). 

(3) $ conserves energy for any H , (q,p) and At for which it is defined. 

(4) $ satisfies Condition 7 

(5) For sufficiently small At, §(H, (q,p), At) is defined, continuous, and con- 
serves volume for (q,p) G U . 

Then for sufficiently small At, there is a function V : K. — > R such 
that if H(q,p) = p\ + V(qi), then <&(H,-,At) is not simultaneously defined, 
continuous, and volume-conserving on U. For each such At, the V constructed 
can be replaced by XV for AG (0, 1) and the same result holds. 

Proof: We will use the hypothesized integrator on M? n to construct an inte- 
grator on M? satisfying the conditions of Theorem 5. Let H : M 2 — > R be 
a given Hamiltonian function. Define H* : M? n — > R by H*(q,p) = H(qi,pi). 
We define the integrator by 

4>{H, ( qi , Pl ), At) = nMH*, ((qi, 0, . . . , 0), ( Pl , 0, . . . , 0)), At). 

Let U C M 2 be given by 

Pi)|((?i,0,...,0),( Pl ,0,...,0))GC/}. 

It is straightforward to check that satisfies the conditions of Theorem 5 on 
U. Thus, by the theorem, we have an arbitrarily small C{f function V such 
that is not simultaneously defined, continuous, and volume-conserving on U 
for H(q 1 ,p 1 )=pj + V(q 1 ). 

Let H*(q,p) =pj/2 + V(q 1 ) for (q,p) G R 2n . Now suppose that $(&*,-, At) is 
defined, continuous, and volume preserving on U . We will derive a contradic- 
tion by showing this implies that (f>(H,-,At) is, in fact, defined, continuous, 
and volume-conserving on U . 

First note that $(if*, •, At) being defined and continuous on U implies that 
(f)(H,-,At) is defined and continuous on U. To check volume conservation, 



10 



note that the Jacobian of the map *, ■, At) has structure 

Ju J12 
/ 

where we have put the variables in order (q 1 ,p 1 , . . . , q n ,p n ) and J n is a 2-by-2 
matrix. Since the determinant of this matrix is 1 by volume-conservation, the 
determinant of Jn must be 1. But Jn is the Jacobian of (j)(H, •, At), so this 
latter map must be area preserving. This contradicts our earlier assumption. 

o 

The second condition states that if the Hamiltonian system consists of n identi- 
cal uncoupled one-degree-of-freedom systems, then the integrator itself should 
consist of n identical uncoupled maps on the state-space of each subsystem. 

Condition 9 If H(q,p) = X^=i M'fr'Pi) then there is an integrator on R 2 
such that 

7T;$(#, (q, P ),At) = <j>(h, (qi,Pi),At) 

for i — 1, ... ,n, for any (q,p) and At for which the integrator is defined. 

Though this is certainly a nice property for the integrator to have (since the 
flow map has the same property) there are many integrators for which it 
does not hold. For example, step-and-project methods may not satisfy this 
condition, even if the underlying one-step method does. 

Theorem 10 Let $ be an integrator for Hamiltonian systems in R 2n . Let 
H(q,p) = HiVi/^- Let U C R 2ri be an open neighbourhood of the set 

{(( qi ,..., qi ),(l,...,l))\ qi eR}. 
Suppose 

(1) $ is consistent (Definition 3). 

(2) $ is computationally reasonable (Definition 2). 

(3) $ conserves energy for any H , (q,p) and At for which it is defined. 

(4) $ satisfies Condition 9 

(5) For sufficiently small At, (q,p), At) is defined, continuous, and con- 
serves volume for (q,p) G U . 

Then for sufficiently small At, there is a function V : R — > R such that 
if H(q,p) = YH=iij>i + y(<Zi)) then $(H,-,At) is not simultaneously defined, 
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continuous, and volume-conserving on U . For each such At, the V constructed 
can be replaced by XV for A G (0, 1) and the same result holds. 



Proof: This theorem is proven analogously to the previous theorem. For any 
H : R 2 -> R we define H* : R 2n — > R by H*(q,p) = EiH(q i:Pi ). We define 
the integrator by 

0(tf, ( qi , Pl ), At) = nMH*, ((gi, . . . , gi), (pi, . . • At). 



We define U G R 2 by 

^ := ((5i, •••,?i),(Pi, •••,?!)) e U}. 



As in the proof of the previous theorem, <fi and U satisfy the conditions of 
Theorem 5. Thus, by the theorem, we have an arbitrarily small function 
V such that <fi is not defined continuous and volume-conserving on U for 
H = pj + V( qi ). 

Let H*(q,p) = Ei(pf/ 2 +^(?i)) for e R2n - Now suppose that $(/?*, •, At) 
is defined, continuous, and volume preserving on U . We will derive a contra- 
diction. 

Now 4>(H, •, At) is defined and continuous on U. To check volume conservation, 
note that the Jacobian of the map Q(H*, •, At) in this case has structure 

Jn ... 
J 22 ■ ■ ■ 
: : "-. 



where we have put the variables in order (qi,pi, ■ ■ ■ , q n ,Pn) and each Ju is 2- 
by-2. Since the determinant of this matrix must be 1 by volume-conservation 
and the determinants of the Ja are identical, the determinant of Jn must be 
±1. As in the proof of the previous theorem this implies (f>(H, •, At) is area 
conserving on U which is a contradiction, o 
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